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Abstract 

Background: As an arborescent and perennial plant, Moso bamboo {Phyllostachys edulis (Carriere) J. Houzeau, synonym 
Phyllostachys heterocycia Carriere) is characterized by its infrequent sexual reproduction with flowering intervals ranging 
from several to more than a hundred years. However, little bamboo genomic research has been conducted on this due to a 
variety of reasons. Here, for the first time, we investigated the transcriptome of developing flowers in IVloso bamboo by 
using high-throughput lllumina GAII sequencing and mapping short reads to the IVloso bamboo genome and reference 
genes. We performed RNA-seq analysis on four important stages of flower development, and obtained extensive gene and 
transcript abundance data for the floral transcriptome of this key bamboo species. 

Results: We constructed a cDNA library using equal amounts of RNA from Moso bamboo leaf samples from non-flowering 
plants (CK) and mixed flower samples (F) of four flower development stages. We generated more than 67 million reads from 
each of the CK and F samples. About 70% of the reads could be uniquely mapped to the Moso bamboo genome and the 
reference genes. Genes detected at each stage were categorized to putative functional categories based on their expression 
patterns. The analysis of RNA-seq data of bamboo flowering tissues at different developmental stages reveals key gene 
expression properties during the flower development of bamboo. 

Conc/usion:\Ne showed that a combination of transcriptome sequencing and RNA-seq analysis was a powerful approach to 
identifying candidate genes related to floral transition and flower development in bamboo species. The results give a better 
insight into the mechanisms of Moso bamboo flowering and ageing. This transcriptomic data also provides an important 
gene resource for improving breeding for Moso bamboo. 
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Introduction 

Transcriptome sequencing is a convenient way of rapidly 
obtaining information on the expressed fraction of a genome [1]. 
Next-generation sequencing (NGS) technology, such as the 
lllumina Solexa, Roche 454 and ABI SOLID platforms, has 
emerged as a cost-effective approach for high-throughput 
sequencing. It has revolutionized biological research by rapidly 
providing genomic and transcriptomic data [2-4] . RNA-Seq refers 
to the whole-transcriptome shotgun sequencing of fragmented 
mRNA or cDNA, resulting in overlapping short reads covering the 
entire transcriptome [5,6]. RNA-Seq analysis is powerful, enabling 



large-scale functional assignments of genes via the assembly of 
large sequenced transcriptome library. With this technique, 
quantitative analysis of gene expression can be easily performed 
without potential bias, allowing a more sensitive and accurate 
profiling of the transcriptome [7,8]. Despite its obvious potentials, 
this combined approach has not yet been applied to Moso bamboo 
{Phyllostachys edulis (Carriere) J. Houzeau, synonym Phyllostachys 
heterocycia Carriere) [9] research. 

Bamboo is one of the most important non-timber forest 
products in the world. About 2.5 billion people depend on 
bamboo, ranging from food to raw materials for construction and 
manufacture [10-12]. Bamboos (Bambusoideae) belong to the 
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monophyletic BEP clade (Bambusoideae, Ehrhartoideae, Pooi- 
deae) in the grass family (Poaceae), and consist of woody and 
herbaceous species. Woody bamboos are arborescent, perennial 
plants characterized by their woody stems and a prolonged 
vegetative phase lasting decades or even longer before flowering 
[13]. In addition, they often flower synchronously and die 
collectively after flowering. Two types of bamboo are usually 
recognized - sympodial bamboos, which form clumps of culms, 
and monpodial bamboos which form groves. Under normal 
circumstances, flowering of monopodial bamboo rarely occurs. 
However, collective death in some bamboo species has occasion- 
ally occurred after flowering over a large area. The reasons for this 
remain unclear, but it is a devastating blow to bamboo resource 
base. Therefore, it is important to determine the specific pathways 
and genes involved in bamboo flowering and flower development. 

Several putative flowering-related genes have been identified 
from certain bamboo species [14—18], and environmental and 
chemical manipulations have been found to induce bamboo 
flowering in vito [19]. Moreover, Zhang et al. studied the 
transcriptomes of developing flowers in Dendrocalamus latifloms 
using lUumina sequencing [20]. 

As a large woody bamboo with high environmental, economic 
and cultural value in Asia [21], Moso bamboo is a perennial plant 
characterized by its infrequent sexual reproduction with flowering 
intervals ranging from several to more than a hundred years. 
Flowering depends on many factors, including environment, 
nutrition, climate and the plants' physiological status. Several 
hypotheses have been proposc'd for the flowering mechanisms 
[22,23]. However, there have been few molecular studies on 
flowering of Moso for various reasons including the difficulty of 
collecting tissue samples, lack of genomic knowledge, infrequent 
sexual reproduction and long periods of time between flowering. 
With the announcement of the genome sequence of Moso bamboo 
[21], it is now possible to identify and determine the molecular 
regulation mechanisms of all functional elements in Moso 
bamboo. The transcriptome represents a comprehensive set of 
transcribed regions throughout the genome. Therefore, under- 
standing the transcriptome dynamics is essential for revealing 
functional elements of the genome and interpreting phenotypic 
variation produced by combinations of genotypic and environ- 
mental factors. 

In the present study, we aimed to investigate the transcriptome 
of developing flowers in Moso bamboo through high-throughput 
lUumina GAII sequencing and mapping short reads to the Moso 
bamboo genome and reference genes. Here, compared with the 
previous study reference [as 21], for the first time we collected 
from the same plant both leaves in the vegetative stage and 
flowering spikelet samples in diflerent floral developmental stages. 
As transcriptome results are affected by the impact of external 
environmental conditions and growth differences among plants, 
transcriptome results from the same plant are more accurate to 
elucidate the molecular mechanism of regulation of flower 
development of bamboo. 

Transcripts from four flower growth phases of the same plant 
were isolated, quantified and sequenced. The transcriptome 
sequences were then blasted against public databases. Subse- 
quently, the annotated genes were clustered into putative 
functional categories using the Gene Ontology (GO) framework 
and grouped into pathways using the Kyoto Encyclopedia of 
Genes and Genomes (KEGG). Finally, we assigned genes putative 
homologs in model species, including Arabidopsis ihalima, Oryza 
sativa, Brachjpodium distcwhjon and other relatives in the grass family 
to determine whether a hierarchy of gene regulation may persist 
between these species. Taken together, we, for the first time. 



comprehensively characterized the molecular basis of the physi- 
ological processes during the flower development of Moso bamboo 
using large-scale high-throughput sequencing, and our results may 
lay a foundation for further gene expression and functional 

genomic studies in bamboos. 

Results and Discussion 

Sequence data quality control and filtering 

In order to get an overview of the Moso bamboo transcriptome, 
we generated cDNA libraries from leaves of non-flowering plants 
(abbreviated as CK) and equally mixed RNA extracted from 
panicles at four different flowering developmental stages (abbre- 
viated as F) using paired-end sequencing by the lUumina platform. 
In total, we acquired more than 67 million reads from each of the 
CK and F samples. Through mapping all these short reads to the 
reference genome of Moso bamboo, we found that about 70% of 
the reads could be unirjuely mapped to the genome. Gene 
coverage was determined as the percentage of a gene covered by 
reads. 15,497 and 167,767 reads were detected from each of the F 
and CK samples respectively, with gene coverage from 90% to 
100% (Figure SI). 

The reads mapping to the reference genome dataset 

To identify the genes corresponding to these clean reads in each 
fibrary, the clean reads were mapped to the reference genes 
expressed in the Moso bamboo genome. Mapping results showed 
that 39,026.185 and 39,386,902 reads from each librar)' were 
perfectly matched to the reference genome while about 22,824,394 
and 23,064,627 reads were perfectly matched to the reference 
genes (Table 1). The reads of unique matched to genome were 
47,847,557 and 47,709,687, and those matched to reference genes 
were 30,640,507 and 30,868,269 in the two libraries. Altogether, 
there were 52,489,347 and 53,035,796 reads matched to the 
reference genome, and 31,843,985 and 32,346,195 reads matched 
to the reference genes. However, as a result of the significant 
sequencing depth of lUumina technology and incomplete annota- 
tion of the Moso bamboo genome, 15,815,789 and 14,582,070 
unmatched reads to the reference genome and 36,461,151 and 
35,271,671 unmatched reads to the reference genes in each library 
were observed. The alignment statistics results of these reads from 
library CK and F are also shown in Table 1. 

GO annotation and KEGG pathway mapping 

Based on sequence homolog)-, a total of 18,309 genes were 
categorized into the three main GO categories (biological process, 
ceUular component and molecular function). The (X) analysis 
showed that the identified genes were involved in flower 
development. In the GO classification, "pollen wall" (6,420, 
94.3%), "exine" (6,420, 94.3%) and "sexine" (5,027, 73.8%) terms 
were dominant, respectively. We also noticed that genes of 
"spindle microtubule" (5,590, 82.1%), "ornithine transport" 
(4,123, 71.5%), "ketone biosynthetic process" (4,119, 71.4%), 
"actin filament severing" (3,087, 53.5%) and "regulation of 
meiotic ceU cycle" (3,080, 53.4%) were weU represented. 
Moreover, additional table SI shows that 1,794 genes were 
annotated as "movement in host emironment" t:ategory, suggest- 
ing a potential connection between environment and flowering. 
Functional classification and pathway assignment were performed 
by die KEGG [24]. A total of 23,522 genes were mapped into 302 
KEGG pathways, representing compound biosynthesis, utiliza- 
tion, metabolism and degradation. The largest category was 
metabolic pathways which included 7,316 genes. Additional 
table S2 shows that transcripts related to genetic information 
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Table 1. Alignment statistics of sample CK and F. 







Tissue 


F 


CK 


Total Reads 


68305136 


67617866 


Total BasePairs 


6147462240 


6085607940 


Map to Genome 


Total mapped reads 


52489347 


53035796 


Perfect match 


39026185 


39386902 


< = 5bp mismatch 


13463162 


1 3648894 


Unique match 


47847557 


47709687 


IVlulti-position match 


4641790 


5326109 


Total unmapped reads 


15815789 


14582070 


Map to Gene 


Total mapped reads 


31843985 


32346195 


Perfect match 


22824394 


23064627 


<= 5bp mismatch 


9019591 


9281568 


Unique match 


30640507 


30868269 


IVlulti-position match 


1203478 


1477926 


Total unmapped reads 


36461151 


35271671 



doi:10.1371/journal.pone.0098910.t001 



processing and environmental information processing were the 
most abundant such as RNA transport (Ko030I3, 3,921 genes, 
16.67%), mRNA surveillance pathway (Ko03015, 3,301 genes, 
14.03%), microbial metabolism in diverse environments 
(K0OII2O, 1,621 genes, 6.89%), plant-pathogen interaction 
(Ko04626, 1,582 genes, 6.73%) and plant hormone signal 
transduction (Ko04075, 1,293 genes, 5.5%). 

Putative Moso bamboo floral development transcription 
factors 

One unique characteristic of Moso bamboo is the switch to 
flowering after a very long period of vegetative growth and the 
undertermined flowering time. In order to compare the gene 
expression patterns between flowering and vegetative tissues, we 
collected panicles and vegetative leaves from non-flowering Moso 
bamboo plants for Illumina sequencing analysis. More than 714 
floral development-related genes were up-regulated in the panicle 
tissues (with at least a 2-fold difference in the expression level in 
panicles compared with the level in leaves), including 476 
flowering genes in the category of transcription factor genes 
(Table 2, Table S3 and Table S4). 

Transcriptional regulation is mediated through the interplay 
between transcription factors and specific-regulatory regions of the 
genome, leading to gene activation or other changes [25]. 
Therefore, some transcription factors may be considered as 
master regulators, since their actions initiate a branching series 
of downstream effects, including the activation of other transcrip- 
tion factors, ultimately resulting in broad changes in the organism 
(as [19]). Our study focused on these key regulatory genes in the 
flower development of Moso bamboo. A total of 476 putative 
transcription factor families were identified in panicle tissues of 
Moso bamboo, including WRKY (110), MADS-Box (38), bHLH 
(34), Dof (15), NAC (61), HSF (4), HMG (4), bZIP (6), HSP20 (28), 
HSP70 (1), MYB (143) and AP2/EREBP (32). 

MADS-box transcription factors contain a DNA binding 
domain conserved among eukaryotes. Many MADS family 
members have been intensely studied in model plants, with the 



floral-organ identity and development proteins falling into the 
MIKC clade [26-28], such as homologs of MADS 1, MADS2, 
MADS5, MADS7, MADS12, MADS14, MADS15 and MADS56 
in Brachypodium distachyon, Oryza sativa Japonica Group and Fargesia 
nitida. MADS family members have been shown to orchestrate 
floral organ specification and development [20,26]. Loss-of- 
function mutants of MADS-box genes have caused changes in 
organ identity. In this study, a total of 38 MADS-Box transcription 
factors were identified in Moso bamboo flowers, which are 
homologs of MADSl, MADS2, MADS4, MADS7, MADSll, 
MADS14, MADS15, MADS17, MADS27, MADS31 and 
MADS58 from ^ea mays, Fargesia nitida, Oryza sativa Indica Group 
and Brachypodium distachyon. Genes OsMADS14, OsMADSlS (rice 
FULl-like gene) is expressed only in inflorescence and developing 
caryopses [29]. OsMADSlS (rice FUL3-like gene) causes an early 
flowering phenotype and early initiation of axUliary shoot 
meristems. It indicates that OsMADSlS acts as a promoter in the 
difierentiation activity of the vegetative shoot [30]. Dreni et al 
indicated that OsMADS3 and OsMADSSS might have a redun- 
dantly function in specifying floral organ identity, including the 
carpel [31]. 

Dof (DNA binding with one zinc finger) transcription factors are 
a group of plant-specific transcription factors containing a single 
Cys2/Cys2-type zinc-finger-like Dof motif [32,33]. There are 37 
Dof transcription factors in the Arabidopsis genome, and all of 
them have a highly conserved DNA-binding domain and 
divergent sequences beyond this domain [33-36]. They are 
involved in diverse functions, such as seed dormancy [37], carbon 
metabolism [38], phytochrome signaling pathway [39], phenyl- 
propanoid metabolism [40] and flowering response [4 1 ,42] . Jci)o/3 
is a circadian clock regulated gene, and it might be involved in the 
regulation of flowering time in Jatropha curcas [43]. Li et al. 
suggested that OsDofl 2 over-expression induces an early flowering 
under long-day (LD) conditions, the transcription levels of Hd3a 
and OsMADS14 were up-regulated [44]. We detected 28 putative 
Dof transcription factors in panicles, which are homologs of Dof3, 
Dof4, DofS, Dofl2 and CDF (Cychng DOF Factor) family, 
revealing similar roles of Dof 3 and Dof 1 2 genes during the flower 
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Table 2. The genes related to transcription factor families. 




Transcription Factor Families 


Gene Numbers 


MADS-Box gene (MADS) 


38 


heat shock transcription factor 20 (HSP20) 


28 


basic helix-loop-helix factor (bHLH) 


34 


basic leucine zipper (bZIP) 


6 


NAM-ATAF-CUC family (NAC) 


61 


heat shock transcription factor 70 (HSP70) 1 


WRKY transcription factor (WRKY) 


110 


MYB transcription factor (MYB) 


143 


DNA-binding with-one finger (Dof) 


15 


AP2/EREBP transcription factor gene (AP2/EREBP) 


32 


Heat Shock Transcription Factor (HSF) 


28 


high-mobility group box (HMG) 


4 


doi:l 0.1 371/journal.pone.009891 0.t002 



development of Moso bamboo. We also identified Hd3a, a master 
floral developmental regulator in rice (with apparent paralogous 
counterparts in Arabidopsis called FT). Peng et al. (2013) suggested 
that drought or other environmental stresses are fianctional in 
regulating MADS14s in the flowering stage of bamboo (as [21]). 
Taken together, an active pathway of Dof-Hd3a-MADS-flowering 
may play an important role during the Moso bamboo flowering. 

WRKY transcription factors are associated with senescence 
[45,46] and stress responses [47]. Members of this family have 
been identified as important downstream components of MAPK 
signaling pathways (as [19]). The involvement of WRXY factors in 
plant defense is also well documented: some of them have been 
shown to confer disease resistance [48,49], triggering the 
expression of defense-related genes in systemic acquired resistance 
[50-53]. They may be induced by signaling hormones, such as 
salicylic acid [54,55], jasmonic acid [56] and gibbereUic acid [57]. 
In the present study, putative homologs of WRKY transcription 
factors exhibited a similar over-representation, with 110 genes 
showing high abundance in panicles of Moso bamboo. 

The basic helix-loop-helix family (bHLH) consists of genes 
regulating various processes of flower development, such as 
controlling the development of carpel margins, as well as the 
morphogenesis of sepals, petals, stamens and anthers in Arabidopsis 
thaliana and Eschscholzia californica [58] . In our study, 34 bHLH-like 
genes were highly expressed in panicle tissues. 

MYB transcription factors contain DNA binding domains, and 
some have been identified as floral developmental regulators [59]. 
In Arabidopsis, MYB21, MYB24 and MYB57 are crucial factors 
for the development of stamen filament. It has been reported that 
R2R3-MYB genes regulate various metabolic pathways, contrib- 
uting to the control of flavonoid biosynthesis, tryptophan 
biosynthesis and participating in phenylpropanoid metabohsm 
[60,61]. In plants, basic leucine zipper motif (bZIP) transcription 
factors regulate a series of processes, including pathogen defense, 
light and stress signaling, seed maturation and flower develop- 
ment. We found that these upstream regulatory stress-responsive 
genes were highly expressed during flowering, suggesting a 
potential connection between drought, pathogen or senescence 
and bamboo flowering. 



Detection of putative genes related to flowering time 
control and flower development 

Previous studies of model plant species showed that the timing 
of floral induction is controlled by sophisticated regulatory 
networks that monitor changes in the environment, including 
autonomous pathway, photoperiod and circadian clock pathway, 
gibberellins pathway, ambient temperature pathway and age 
pathway [62,63]. The unusual and infrequent nature of bamboo 
flowering has attracted the curiosity of scientists and laypeople for 
centuries [64] . However, little research has been conducted about 
the molecular mechanism of bamboo flowering due to the 
difficulty of collecting flowering samples. Through the comparison 
of Moso bamboo genes found in this study with the NCBI and 
Uniprot databases, we identified at least 238 genes as homologs of 
known flowering-related genes from other plants (Table S3). 
These genes were compared with flowering genes from Oryza 
saliva, Hordeum vulgare, Brachypodium distachyon, Sorghum bicolor and 
Arabidopsis thaliana. However, we found that the genes employed in 
typical flowering promotion pathways (such as those autonomous 
pathway, ambient-temperature) and floral pathway integrator 
(FPI) genes (as [63,65]) were not highly expressed in these floral 
tissues in bamboo. Only late elongated hypocotyl [LHT), phyto- 
chrome-interacting facor 3 {PIF3), crytochrome- 1 [CRYl), GI- 
GANTEA {GI) and CONSTANS (CO) genes were detected in 
Moso bamboo for the photoperiod and circadian clock pathway. 
We did not detect any genes homologous to components of the 
circadian clock and photoperiod pathway, such as early flowering 
3 [ELF3), early flowering 4 {ELF4), circadian clock associated 1 
(CCAl), timing of CABl {TOCl), pseudo response regulator {PRR5, 
PRR7 and PRR9) and phtyochrome {PHTA, PHTB) [66]. The CO 
gene plays a key role in integrating light and temporal information, 
and it is essential for determining the flowering time [67]. The CO 
protein consists of two zinc fingers, including a N-terminal B- 
domain mediating protein-protein interactions and a C-terminal 
CCT domain for nuclear localization [68]. In Arabidopsis, CO, as 
a transcription factor, promotes flowering by inducing the 
expression of flowering locus T (FT) [69]. Rice orthologues of 
the Arabidopsis genes CO and FT haye been identified as Hdl and 
Hd3a, respectively; and Hdl controls rice flowering by regulating 
Hd3a (as [57]). In the present study, we found that 41 putative CO 
transcription factors were involved in flower development. 
However, most of them were down-regulated in floral tissues of 
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Moso bamboo. The upstream floral repressor factors, such as 
LHY and PIF3, were highly expressed, whereas the FPI genes (G/ 
and CRTl) were down-regulated, which might result in the low 
expression of CO genes. Low expression of CO genes suggested that 
the flower initiation was more independent of these traditional 
promotion pathways in bamboo flowering. 

Genes involved in gibberellins pathway have high expression, 
which is known as the promotion pathways in sophisticated 
regulatory networks. GA regulates diverse developmental process- 
es during plant growth, ranging from seed germination, leaf 
expansion, and stem elongation to floral initiation [70]. Putative 
homologs of GA-signaling pathway genes were highly expressed, 
including the gibbereUin response modulator (GAMYB) [71] and 
gibberellin receptors [GIDl, GID2) [72]. GibbereUin triggers 
degradation of the DELLAs (such as GAI, RGA, and RGL78-80) 
and up-regulates MYB21, MYB24 and MYB57, which are essential 
for stamen filament growth (as [59]). High expression of these 
genes indicated the connection between GA and bamboo 
flowering. 

Unlike the flowering pathway genes, over 252 genes involved in 
the response to drought stress or to other correlative stresses (such 
as oxidative stress) were highly expressed. Some FMI-related genes 
and their upstream regulatory drought-responsive genes were 
highly expressed during flowering, such as DoJ3, Dqf5, DofJ, Dofl 0, 
Dofl2 and bZIP transcription factors. Putative homologs of FMI 
genes [73-75], such as API and AP2, MADSl, MADS14 and 
MADS15, were also highly expressed. High expression of these 
above-mentioned genes indicated a potential connection between 
drought or other environmental stress and bamboo flowering. 

Changes in gene expression profiles during different 
flower developmental stages 

DGEs (Differentially Expressed Genes) generates absolute 
rather than relative gene expression measurements and avoids 
many inherent limitations of microarrays (as [8]). DGEs were used 
to analyze the gene expression during four stages of flower 
development in Moso bamboo. Four cDNA libraries (Fl, F2, F3 
and F4) were sequenced, resulting in between 1 0 to 11 million 
high-quality reads. The read sequences were mapped to Moso 
bamboo reference transcriptome database. About 8,412,664 to 
9,367,437 reads were mapped to a gene in the reference genome 
database, and up to 80% of the sequences in transcriptome could 
be unequivocally identified to the Moso bamboo genome. About 
5,185,207 to 5,820,078 reads were mapped to a gene in the 
reference gene database, and up to 79% of the sequences in 
transcriptome could be perfect matched to the Moso bamboo 
genes (Table 3). Additional Figure S2 shows the quality assessment 
of reads, sequencing saturation analysis, randomness assessment 
and gene coverage. We identified differentially expressed genes 
between samples using a previously developed algorithm by Audic 
and Claverie [76] . As expected, the majority of gene expression 
changes occurred between the CK and F1/F2/F3/F4 samples, 
with more up-regulated genes observed (Figure 1 & Table S5, S6, 
S7, and S8). Moreover, about 1% (Table S9) of the up-regulated 
genes were found to be orphan sequences, which were not 
annotated to the Gene Ontology database (http://www. Geneon- 
tology.org/) with nr annotation, but the expression level of genes 
were significantly different, so we speculate that these genes may 
take part in other unique processes and pathways involving in the 
flower development of Moso bamboo. 



up-regulated "down-regulated 




CK CK CK 

vs. vs. vs 

F1 F2 F3 



Figure 1. Changes in gene expression profile among tKie 
different flower developmental stages. The number of up- 
regulated and down-regulated genes between CK and Fl; CK and F2; 
CK and F3; CK and F4; Fl and F2; Fl and F3; Fl and F4; F2 and F3; F2 
and F4; F3 and F4 are summarized. DEGs: Differentially Expressed 
Genes. 

doi:1 0.1 371/journal.pone.009891 0.gOOl 

Expression levels of putative genes involved in flower 
development 

A total of 80 genes were annotated as flowering regulating 
pathway in Moso bamboo transcriptome (Table SIO). To find 
whether these genes were putatively involved in flowering 
regulating pathway, expression profiles of these genes were 
compared with CK sample by DGEs (shown in Figure 2). From 
CK to Fl and from CK to F2, MADS, Dof, GIDl, GID2, LHT and 
MTB genes were up-regulated. Some of these genes exhibited a 
decreased expression over time. The expression oi Dof genes was 
significandy decreased in the anthesis and embryo formation 
stages, suggesting that Dof genes played a greater role in the early 
stages of flower development. Flowering inhibitors, such as LHT, 
PIF3 genes, were up-regulated during flower developmental 
periods, which might result in the down-regulated expression of 
CO genes in floral tissues. Low expression of CO genes and high 
expression of stress-response genes and FMI genes suggested that 
the activation of FMI was more independent of the traditional 
flowering regulating pathway in Moso bamboo flowering. 

Gene validation and expression analysis 

To confirm experimentally that the genes obtained from clean 
reads and computational analysis were indeed expressed, 8 
putative genes related to flower development were chosen for 
RT-PCR and qRT-PCR analysis (Figure 3 and Table Sll). 

In the RT-PCR analysis, every selected gene was PCR positive 
with a single band at the calculated size (data not shown). 
According to the qRT-PCR results (Figure 3 and table 4), CO gene 
were down-regulated in floral tissues of Moso bamboo, whereas 
LHT and PIF3, which are the upstream floral repressor factors, up- 
regulated during the flower development. Drought-responsive 
genes such as Dof were highly expressed in floral tissues, especially 
in the early stage of flowering; GIDl expression levels were similar 
to Dof in floral tissues. The expression levels oi MTB, MADS- 14, 
AdADS-15 in floral tissues were higher than those in non-flowering 
leaves. The results of qRT-PCR expression analysis matched the 
putative functions of lUumina analysis. 

Conclusions 

Littie genetic or genomic research has been performed into the 
fascinating mystery of bamboo flowering despite the economic 
importance of several bamboo species (as [20]). We investigated 
the sequences and transcript abundance levels at different 
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Table 3. Statistics of RNA-seq sequencing. 





Tissue 


F1 


F2 


F3 


F4 


Total Reads 


10424052 


1 1 1 27400 


11046632 


10828810 


Total BasePairs 


510778548 


545242600 


541284968 


530611690 


Map to Genome 


Total mapped reads 


8768572 


9367437 


9101401 


8412664 


Perfect match 


7184271 


7662051 


7467544 


6787731 


<= 2bp mismatch 


1584301 


1 705422 


1633857 


1 624933 


Unique match 


7458759 


7942000 


7759666 


7199619 


Multi-position match 


1309813 


1425437 


1341735 


1213045 


Total unmapped reads 


1655480 


1 759963 


1945231 


2416146 


Map to Gene 


Total mapped reads 


5626520 


5820078 


5334090 


5185207 


Perfect match 


4452869 


4604128 


4226029 


4028547 


< = 2bp mismatch 


1173651 


1215950 


1108061 


1 1 56660 


Unique match 


5030921 


5208332 


4809727 


4680929 


Multi-position match 


595599 


611746 


524363 


504278 


Total unmapped reads 


4797532 


5307322 


5712542 


5643603 



doi:10.1371/journal.pone.0098910.t003 



developmental stages of Moso bamboo flowers using the lUumina 
GAII platform. We showed the feasibility of using combination of 
transcriptome sequencing and RNA-Seq analysis to study the 



genes involved in floral development in Moso bamboo. More than 
714 floral development-related bamboo genes were highly 
expressed in the panicle tissues. Low expression of CO and FPI 
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Figure 2. Clustering of genes involved in bamboo flowering expression profiles at four different flower developmental stages. 

doi:1 0.1 371/journal.pone.009891 0.g002 
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Figure 3. The expression profiles of 8 selected genes from flowering tissues in different flower developmental stages and leaves of 
non-flowering of Moso bamboo. The transcript levels were normalized to that of TIP41 (tonoplast intrinsic protein 41), and the level of each gene 
in the control was set at 1.0. Error bars represent the SD for three independent experiments. CO (CONSTANS), GIDl (gibberellin receptor GID1), DOF 
(Dof zinc finger protein 12), MYB (MYB transcription factor), FIP3 (phytochrome- interacting factor 3), MADS-14 (MADS-box gene-14), MADS-15 (IVIADS- 
box gene-15), LHY (late elongated hypocotyl). 
doi:1 0.1 371 /journal.pone.009891 0.g003 



genes suggested that FMI was more independent of these known 
promotion pathways in Moso bamboo flowering. Moreover, high 
expression of genes involved in stress-responsive pathways and 
GA-signaling pathway indicated the potential connection between 
adversity stress, GA and bamboo flowering (shown in Figure 4). 
The dataset will provide an important resource foundation for 
future genetic or genomes studies on bamboo species, and will help 
to give better insight into the mechanism of the flower 
development of Moso bamboo. 

Materials and Methods 

Ethics Statement 

All plant protocols were reviewed and approved by the 
International Center for Bamboo and Rattan and Key Laboratory 
of Bamboo and Rattan Science and Technology, State Forestry 



Administration. All necessary permits were obtained for the field 
studies from Guangxi Provincial Academy of Forestry and GuUin 
Forestry Bureau, Dajing County in Guangxi Provence. The field 
work conducted for sampling did not aflect the local ecology and 
did not involve endangered or protected species. 

Sample preparation and RNA extraction 

Moso bamboo samples of non-flowering plants and flowering 
plants at different flowering developmental stages (floral bud 
formation, inflorescence development, anthesis and embryo 
formation stages) were collected in Dajing County, GuUin (E 
110°17'-110°47'; N 25°04'-25°48') in Guangxi Zhuang Autono- 
mous Region from April to August, 2012 (see Figure 5 for details). 
Four developmental stages were defined, based on the anatomical 
structure of floral organs: Fl, F2, F3, and F4. At the first stage (Fl), 
when the floral bud had begun to form, the plant transits from the 
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Figure 4. A hypotKiesized pathways in regulation of flowering 
in iVIoso bamboo. Yellow and red oval Indicate that the involved 
genes are more highly expressed in the floral tissues, whereas blue oval 
indicate that the genes are down-expressed and green oval indicate 
that the genes are not activated. Purple dashed arrows represent 
pathways that were not used during flowering. Black arrow represents 
promotion. Green bar represents suppression. 
doi:1 0.1 371/journal.pone.009891 0.g004 

vegetative to the reproductive stage. At the second stage (F2), the 
inflorescence axis continued to extend, lateral buds started to 
diflerentiate, and panicles grow from the plant but without 
flowering. At the third stage (F3), anthesis, flowers with both pistils 
and stamens emerged from glumes; flowering spikelets were 
collected at this stage. At the last stage (F4), pistils and stamens had 
withered, and the embryo was forming. 

Samples were immediately frozen in liquid nitrogen and stored 
at — 80°C until further analysis. Total RNA was extracted using 
the Trizol reagent (Invitrogen, USA). The quality and purified 
RNA was initially assessed on an agarose gel and NanoDrop 8000 
spectrophotometer (NanoDrop, Thermo Scientific, Germany), 
and then the integrity of RNA samples was further evaluated using 
an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). 

Library preparation for transcriptome analysis 

Samples of two different growth periods were prepared for 
transcriptome analysis using the lUumina's kit (lUumina, San 
Diego, CA, USA) according to the maiwfacturer's instructions. 
Briefly, mRNA was extracted using magnetic beads with Oligo 
(dT) (for eukaryotes) or by removing rRNAs from the total RNA 
(for prokaryotes). Purified mRNA was then mixed with the 
fragmentation buffer and fragmented into short fragments, which 
were used as templates for cDNA synthesis. Purified short cDNA 
fragments were resolved in EB buffer for end reparation and single 
nucleotide A (adenine) addition, and then they were connected 
with adapters. Subsequently, suitable cDNA fragments were 
selected for the PGR amplification as templates. During the QC 
steps, quantification and qualification of the sample library were 
assessed using an Agilent 2100 Bioanaylzer and a ABI StepOne- 
Plus Real-Time PGR System. Finally, the library was sequenced 
using lUumina HiSeq 2000 at Beijing Genomics Institute (BGI) in 
Shenzhen, Ghina. 
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Figure 5. Examples of sampled floral tissues. A1-A2: Representatives of flowers collected to produce the sample for F1 (floral bud formation 
stage). B1-B2: Representatives of flowers collected to produce the sample for F2 (inflorescence growing stage). C1-C2: Representatives of flowers 
collected to produce the sample for F3 (bloom stage, flowers with both pistils and stamens emerging from glumes). D1-D3: Representatives of 
flowers collected to produce the sample for F4 (embryo formation stage). 
doi:1 0.1 371 /journal.pone.009891 0.gOOS 



Sequence data quality control and filtering 

As Moso bamboo genomic resources were available as the 
reference, the original image data was transferred into sequence 
data as raw data or raw reads via base calling after the libraries 
had been sequenced. Then, the clean reads were selected via 
analysis and assessment of base composition and quality, and 
filtering of raw reads by using lUumina PIPILINE software. Since 
the algorithms used in transcriptome construction of the reads 
provided by the lUumina platform may be severely inhibited by 
sequencing errors, a stringent cDNA sequence filtering process was 
employed to select clean reads. Firstly, reads with adaptors were 
removed. Secondly, reads in which unknown bases comprised 
more than 5 % of the total were removed. Finally, low quality reads 
in which the percentage of low quality bases (base quality £20) was 
over 30%, were removed. AH the clean reads were deposited in the 
National Center for Biotechnology Information (NCBI) and can 
be accessed in the Short Read Archive (SRA) linking to BioProject 
accession number: SRRl 187864, SRRl 185317. 

Mapping short reads to the Moso bamboo genome and 
reference gene 

The Moso bamboo genome and reference gene set were 
downloaded from the National Center of Genome Research 
(http://202.127.18.221/bamboo/index.php), and NCBI site 
(http://www.ncbi.nlm.nih.gov/ nuccore/FO203/436, http://www. 
ncbi.nlm.nih.gov/ nuccore/FO203437, http://www.ncbi.nlm.nih. 
gov/ nuccore/FO203443, http:/ / www.ncbi.nlm.nih.gov/ nuccore/ 
FO203444, http://www.ncbi.nlm.nih.gov/nuccore/FO203447, 
http://www.ncbi.nlm.nih.gov/nuccore/FO203448, http://www. 
ncbi.nlm.nih.gov/nuccore/FO203439). 

After reads containing sequencing adapters and low-quality 
reads (reads containing Ns>5) were removed, the remaining reads 
were aligned to the Moso bamboo genome using SOAPaligner/ 
SOAP2 [77], allowing up to two mismatches in each segment 



alignment. For annotation, transcripts were aligned to four public 
databases [no redundant protein database (Nr) in NCBI, non- 
redundant nucleotide database (Nt) in NCBI, Swiss-Prot and 
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway 
database] by BLAST (E-value <10-5). Gene ontology (GO) 
classification was analyzed by the Blast2GO software (v2.5.0) 
based on Nr annotation. 

DGEs library preparation and sequencing 

Tag library preparation for samples of four different flowering 
developmental periods (floral bud formation, inflorescence devel- 
opment, anthesis and embryo formation stages) was performed in 
parallel using the lUumina gene expression sample preparation kit. 
Briefly, mRNA was extracted using magnetic oligo (dT) beads 
from total RNA. Purified mRNA was then mixed with the 
fragmentation buffer and fragmented into short fragments, which 
were used as templates for cDNA synthesis. Purified short cDNA 
fragments were resolved in EB buffer for end reparation and single 
nucleotide A (adenine) addition, and then they were connected 
with adapters. Subsequently, suitable cDNA fragments were 
selected for the PCR amplification as templates. During the Q_C 
steps, quantification and qualification of the sample library were 
assessed using an Agilent 2100 Bioanaylzer and an ABI 
StepOnePlus Real-Time PCR System. Finally, the library was 
sequenced using lUumina HiSeq 2000 at Beijing Genomics 
Institute (BGI) in Shenzhen, China. 

Evaluation of DGEs libraries 

A statistical analysis of the frequency of each tag in different 
cDNA libraries was performed to compare the gene-expression in 
different developmental stages. Statistical comparison was con- 
ducted using the method previously described by Audio and 
Claverie [76]. After the expressional abundances in each library 
were normalized to transcript per mUlioii (RPKM), the most 
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difTerentially regulated genes (differentially expressed genes, 
DEGs) were selected. The significance of differential transcript 
abundance was determined using the FDR (False Discovery Rate) 
control method [78] to justify the p-values, and only genes with an 
absolute fold change of &2 and a FDR significance score of ^ 
0.001 were included in subsequent analysis. In addition, p value 
between two samples was determined using the following formula: 



p{y\x) = 



'N2 

m 



Where Nl and N2 represent the total number of reads mapped 
to genome in each sample, and x and y represent the number of 
reads mapped to a common gene in phase 1 and phase 2, 
respectively. 

Clustering of candidate gene expression profiles 

Hierarchical clustering of log-transformed expression data was 
carried out using the Cluster 3.0 and Treeview programs [79]. 
Correlations between gene clusters were determined using 
Pearson's correlation. Heat maps were constructed using the 
University of Toronto BAR Heatmapper tool (http://www.bar. 
utoronto.ca/ ntools/cgi-bin/ ntools_heatmapper.cgi). 

Quantitative real-tlnne PGR (qRT-PCR) 

To evaluate the validity of lUumina analysis and assess the 
expression profiles in terms of specific mRNA abundances, 8 
putative genes were selected and detected by qRT-PCR, and to 
further investigate the expression profiles of these genes, qRT- 
PCR of flowering tissues at different flowering developmental 
stages were performed separately using leaves of no-flowering 
plants as the control. Reverse transcription reactions were 
performed using 2 mg of RNA by M-MLVRT (Promega, USA) 
according to the manufacturer's instructions. Sequences of 8 
selected genes w(;re obtained from the Moso bamboo genome 
database (http://www.ncgr.ac.cn/bamboo). Primers, picked by 
using the Primer 3 software (http://www.genome.wi. mit.edu/cgi- 
bin/primer/primer3.cgi), as shown in Additional file 13. Tono- 
plast intrinsic protein {TIP41), cited from Fan et al. [80], were used 
as the internal housekeeping gene control. Real-time PGR 
reactions were carried out with LightCycler480 System (Roche, 
USA) using SYBR Premix EX Taq kit (Roche, USA).The 20 jtl 
reaction mixture contained 0.4 ^1 (10 \\M) of each primer and 2 ^ll 
(50 ng) cDNA and 10 |a,l SYBR Green I Master according to the 
manufacturer's instructions. Amplification reactions were per- 
formed as the following: 95°C for 10 s, 60°C for 10 s, and 72°C 
for 20 s. All reactions were performed in triplicate, both technical 
and biological. Data was analyzed using Roche manager software. 

Additional Information 

Accession code 

The lUumina HiSeq 2000 sequencing data for Moso bamboo 
flower transeriptome and RNA-seq have been deposited in NIH 



Short Read Archive database under the accession number 
SRRl 187864 and SRRl 185317. 
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